# ----------------------------------------------------------------------
#  Project: Affective Polarization - Experiment
#  Last updated: Sat Dec 21 18:01:04 2019
#  Purpose: Clean data and generate basic tables and figures
#  Outputs: 
#  Machine: Chagai's macbook pro
# ----------------------------------------------------------------------

# ----------------------------------------------------------------------
# load all relevant themese and packages
# ----------------------------------------------------------------------
rm(list=ls())
library("tidyverse")
library("ggplot2")
library("stargazer")
library("xtable")
library("texreg")
library("reshape2")
library("wesanderson")
library("RItools")
library("estimatr")
library("fastDummies")
library("ggpubr")
library("estimatr")
library("mediation")
library("car")
library("table1")
library("standardize")
library("tidyverse")
library("readstata13")
library("naniar")
library("lfe")

# ----------------------------------------------------------------------
# Read data
# ----------------------------------------------------------------------

main_data <- read_csv("Data for analysis/Survey_analysis.csv")

# ----------------------------------------------------------------------
# Create Figure 2
# ----------------------------------------------------------------------

# Recode Treatment Manipulation check higher # Higher likelihood of unity gov
# This results in a 1-5 scale where 5 is most likely unity
main_data <- main_data %>% 
  mutate(.,
         manipulation_check_1 = 6- (manipulation_check_1-1))

# Create a function to generate mean + ci
calc_range <-function(x){
  n <- length(x)
  mean <-mean(x, na.rm = T)
  high <- sd(x, na.rm = T)*1.96/sqrt(n) + mean
  low <- mean - (sd(x, na.rm = T)*1.96/sqrt(n))
  result <-data.frame(high, mean, low)
  names(result) <-c("high", "mean", "low")
  return(result)
}


unity_cond <- main_data %>% filter(treat == "Unity Government")
narrow_cond <- main_data %>% filter(treat == "Narrow Government")


ttest <- t.test(unity_cond$manipulation_check_1, narrow_cond$manipulation_check_1)


###Likelihood of Unity Government across conditions (manipulation check)
# Subset data
unity_mean_mc <- calc_range(unity_cond$manipulation_check_1) %>% 
  mutate(Condition = "Unity Government")

narrow_mean_mc <- calc_range(narrow_cond$manipulation_check_1) %>% 
  mutate(Condition = "Narrow Government")


t.test(unity_cond$manipulation_check_1, narrow_cond$manipulation_check_1)
# Bind models together
plot_means <- rbind(unity_mean_mc,narrow_mean_mc)

limits_means <- aes(ymax = plot_means$high,
                    ymin = plot_means$low)
scale_means <- c("Unity Government", "Narrow Government")
ggplot(plot_means, aes(x = Condition, y = mean)) +
  geom_errorbar(aes(ymin=low, ymax=high), width = 0.25,
                size = 0.7, alpha = 1, color = "darkblue") +
  geom_point(color = "red", size = 2, alpha = 0.6)+
  labs(x = "Condition",
       y = "Likelihood of Unity Government",
       title = "") +
  annotate("text", x = 2.2, y =3.1, 
           label = "paste(italic(p), \" < 0.05\")", parse = TRUE, 
           family = "Times", size = 3) +
  annotate("text", x = 2.2, y =3.12, 
           label = "paste(italic(Delta), \" = 0.15\")", parse = TRUE, 
           family = "Times", size = 3.1) +
  theme(text = element_text(size = 10, family = "Times"),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(size = 10),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))
ggsave("Figures/manipulation.pdf",  width = 4.85, height = 3.4)


# ----------------------------------------------------------------------
# Create Table 3
# ----------------------------------------------------------------------

# Main model using SD polarization without center

# Polar base
sd_1 <- lm(general_socd_polar_nc ~ unity,
           data = main_data)
summary(sd_1)


# Polar controlling for PM
sd_1_pm <- lm(general_socd_polar_nc ~ unity + pm_eng,
              data = main_data)
summary(sd_1_pm)


# Polar controlling for PM and demographics
sd_1_pm_cov <- lm(general_socd_polar_nc ~ unity + pm_eng +
                    as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                    as.factor(relig) + as.factor(locality) + as.factor(vote),
                  data = main_data)
summary(sd_1_pm_cov)

# In party affect

sd_in <- lm(general_socd_nc_in ~ unity,
            data = main_data)
summary(sd_in)

# Out party affect

sd_out <- lm(general_socd_nc ~ unity,
             data = main_data)
summary(sd_out)

# Create table of main results

stargazer(sd_1, sd_1_pm, sd_1_pm_cov, sd_in, sd_out,
          out= "Tables/experiment_sd.tex",  style="qje",
          title="Effects of Unity Government on Polarization and Party Affect (Social Distance)",
          dep.var.labels = c("Polarization", "In-Party Affect", "Out-Party Affect"),
          keep = "unity",
          #omit = c(3:31),
          covariate.labels = c("Unity"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
          no.space = T,
          omit.table.layout = "n",
          star.cutoffs = NULL,
          font.size = "small",
          label = "tab:exp_sd",
          star.char = c("", "", ""),
          add.lines = list(
            c("", "", "", "", "", "" ),
            c("PM Control", "No", "Yes", "Yes", "No", "No" ),
            c("Demographic Controls",  "No", "No", "Yes", "No", "No"),
            c("Center Voters", "No", "No", "No", "No", "No")))

# ----------------------------------------------------------------------
# Create Table A22
# ----------------------------------------------------------------------

main_data <- main_data %>%
  mutate(.,
         sex_1 = case_when(
           sex == 1 ~ "Male",
           sex == 2 ~ "Female"),
         age_1 = case_when(
           age == 1 ~"18-22",
           age == 2 ~ "23-29",
           age == 3 ~ "30-39",
           age == 4 ~  "40-49",
           age == 5 ~ "50-70"),
         ethnicity_1 = case_when(
           ethnicity == 1 ~ "Ashkenazi",
           ethnicity == 2 ~ "Mizrahi",
           ethnicity == 3 ~ "Russian",
           ethnicity == 4 ~ "Ethiopian",
           ethnicity == 5 ~ "Other"),
         religiosity_1 = case_when(
           relig == 1 ~ "Secular",
           relig == 2 ~ "Traditional",
           relig == 3 ~ "Religious",
           relig == 4 ~ "Haredi"),
         locality_1 = case_when(
           locality == 1 ~ "Jerusalem",
           locality == 2 ~ "Tel Aviv",
           locality == 3 ~ "North",
           locality == 4 ~ "South",
           locality == 5 ~ "Sharon"))

main_data <- main_data %>% 
  mutate(.,
         female = ifelse(sex_1 == "Female", 1, 0),
         male = ifelse(sex_1 == "Male", 1, 0),
         age_a = ifelse(age_1 == "18-22", 1, 0),
         age_b = ifelse(age_1 == "23-29", 1, 0),
         age_c = ifelse(age_1 == "30-39", 1, 0),
         age_d = ifelse(age_1 == "40-49", 1, 0),
         age_e = ifelse(age_1 == "50-70", 1, 0),
         ethnicity_a = ifelse(ethnicity_1 == "Ashkenazi", 1, 0),
         ethnicity_b = ifelse(ethnicity_1 == "Mizrahi", 1, 0),
         ethnicity_c = ifelse(ethnicity_1 == "Russian", 1, 0),
         ethnicity_d = ifelse(ethnicity_1 == "Ethiopian", 1, 0),
         ethnicity_e = ifelse(ethnicity_1 == "Other", 1, 0),
         religiosity_a = ifelse(religiosity_1 == "Secular", 1, 0),
         religiosity_b = ifelse(religiosity_1 == "Traditional", 1, 0),
         religiosity_c = ifelse(religiosity_1 == "Religious", 1, 0),
         religiosity_d = ifelse(religiosity_1 == "Haredi", 1, 0),
         locality_a = ifelse(locality_1 == "Jerusalem", 1,0),
         locality_b = ifelse(locality_1 == "Tel Aviv", 1,0),
         locality_c = ifelse(locality_1 == "North", 1,0),
         locality_d = ifelse(locality_1 == "South", 1,0),
         locality_e = ifelse(locality_1 == "Sharon", 1,0))


disc_table <- dplyr::select(main_data,
                            c(female, male, age_a, age_b, age_c, age_d, age_e,
                              ethnicity_a, ethnicity_b, ethnicity_c, ethnicity_d, ethnicity_e,
                              religiosity_a, religiosity_b, religiosity_c, religiosity_d,
                              locality_a, locality_b, locality_c, locality_d, locality_e, manipulation_check_1,
                              general_therm_polar_nc, general_socd_polar_nc, general_therm_nc,
                              general_therm_nc_in, general_socd_nc, general_socd_nc_in))


stargazer(as.data.frame(disc_table),
          covariate.labels = c("Female",  "Male", "18-22", "23-29", "30-39", "40-49", "50-70", 
                               "Ashkenazi", "Mizrahi","Russian","Ethiopian", 
                               "Other", "Secular", "Traditional", "Religious", 
                               "Haredi", "Jerusalem", "Tel Aviv", "North", "South", "Sharon", "Manipulation",
                               "Polarization (Therm)", "Polarization (Socd)",
                               "Affect out-party (Therm)", "Affect in-party (Therm)",
                               "Affect out-party (Socd)", "Affect in-party (Socd)"),
          title = "Descriptive Statistics - Survey Respondents (Study II)",
          label = "tab:dsc_experiment",
          style = "qje",
          notes = c("Therm refers to feelings thermometer, Socd refers to social disctance scale."),
          notes.append = T,
          notes.align = "l",
          out="Tables/dscrpt_exp.tex")


# ----------------------------------------------------------------------
# Create Figure A4
# ----------------------------------------------------------------------

main_data <- main_data %>% 
  mutate(.,
         treat1 = case_when(
           treat == "Narrow Government" ~ "Narrow \nGovernment",
           treat == "Unity Government" ~ "Unity \nGovernment"
         ))

# Plot assingment to treatment
ggplot(main_data, 
       aes(x = treat)) +
  geom_bar(aes(color = pm_eng,
               fill = pm_eng),
           alpha = 0.7,
           position = "dodge") +
  scale_color_manual(values = c("firebrick2", "dodgerblue2")) +
  scale_fill_manual(values = c("firebrick2", "dodgerblue2")) +
  labs(x = "",
       y = "Frequency",
       color = "",
       fill = "")+
  theme(text = element_text(size = 12, family = "Times"),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(size = 12),
        plot.caption = element_text(size = 12, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))
ggsave("Figures/treatment_n.pdf",  width = 6, height = 6)


# ----------------------------------------------------------------------
# Create Table A23
# ----------------------------------------------------------------------

table1(~ sex_1 + age_1 + ethnicity_1 + locality_1 + religiosity_1 + left_right_1|
  treat, data=main_data, excel=1)


# ----------------------------------------------------------------------
# Create Table A24
# ----------------------------------------------------------------------

# Main model using therm polarization without center
# Polar base
therm_1 <- lm(general_therm_polar_nc ~ unity,
              data = main_data)
summary(therm_1)

# Polar controlling for PM
therm_1_pm <- lm(general_therm_polar_nc ~ unity + pm_eng,
                 data = main_data)
summary(therm_1_pm)


# Polar controlling for PM and demographics
therm_1_pm_cov <- lm(general_therm_polar_nc ~ unity + pm_eng +
                       as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                       as.factor(relig) + as.factor(locality) + as.factor(vote),
                     data = main_data)
summary(therm_1_pm_cov)

# In party affect

therm_in <- lm(general_therm_nc_in ~ unity,
               data = main_data)
summary(therm_in)

# Out party affect

therm_out <- lm(general_therm_nc ~ unity,
                data = main_data)
summary(therm_out)

# Create table of main results

stargazer(therm_1, therm_1_pm, therm_1_pm_cov, therm_in, therm_out,
          out= "Tables/experiment_therm.tex",  style="qje",
          title="Effects of Unity Government on Polarization and Party Affect (Thermometer)",
          dep.var.labels = c("Polarization", "In-Party Affect", "Out-Party Affect"),
          keep = "unity",
          #omit = c(3:31),
          covariate.labels = c("Unity"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
          no.space = T,
          omit.table.layout = "n",
          star.cutoffs = NULL,
          font.size = "small",
          label = "tab:exp_therm",
          star.char = c("", "", ""),
          add.lines = list(
            c("", "", "", "", "", "" ),
            c("PM Control", "No", "Yes", "Yes", "No", "No" ),
            c("Demographic Controls",  "No", "No", "Yes", "No", "No"),
            c("Center Voters", "No", "No", "No", "No", "No")))

# ----------------------------------------------------------------------
# Create Figure A5
# ----------------------------------------------------------------------

# Number of obs for which both SD is positive and therm is negative 14
cor(main_data$general_socd_polar_nc, 
    main_data$general_therm_polar_nc,
    use = "complete.obs")
# Correlation of both measures 60

ggplot(main_data, aes(general_socd_polar_nc, general_therm_polar_nc)) +
  geom_point(alpha = 0.1) +

geom_smooth(method=lm)+
  labs(x = "Social Distance Polarization (No Center)",
       y = "Thermometer Polarization (No Center)",
       caption = expression(rho*"=0.6"))+
  theme(text = element_text(size = 10, family = "Times"),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(size = 10),
        plot.caption = element_text(size = 10, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))
ggsave("Figures/corr_outcomes.pdf", width = 6.5, height = 4)

# ----------------------------------------------------------------------
# Create Table A25
# ----------------------------------------------------------------------

#create variables for likud, blue-white, right, left, and center voters
main_data <- main_data %>% 
  mutate(.,
         likud = ifelse(party_vote == "Likud", 1, 0),
         bw = ifelse(party_vote == "Blue-White", 1, 0),
         left_ideology = ifelse(left_right_1 > 4, 1,0),
         right_ideology = ifelse(left_right_1 < 4, 1,0),
         center_ideology = ifelse(left_right_1 == 4, 1, 0))

#Affect towards Arabs

# unity treatment on affect towards arabs
sd_arabs <- lm(sdis_arb ~ unity,
               data = main_data)
summary(sd_arabs)

sd_arabs_control<- lm(sdis_arb ~ unity+ as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                        as.factor(relig) + as.factor(locality),
                      data = main_data)
summary(sd_arabs_control)

# unity treatment on affect towards arabs conditional on right
sd_arabs_right <- lm(sdis_arb ~ unity*right_ideology+ as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                       as.factor(relig) + as.factor(locality),
                     data = main_data)
summary(sd_arabs_right)

# unity treatment on affect towards arabs conditional on left
sd_arabs_left <- lm(sdis_arb ~ unity*left_ideology+ as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                      as.factor(relig) + as.factor(locality),
                    data = main_data)
summary(sd_arabs_left)

# unity treatment on affect towards arabs conditional on center
sd_arabs_center <- lm(sdis_arb ~ unity*center_ideology+ as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                        as.factor(relig) + as.factor(locality),
                      data = main_data)
summary(sd_arabs_center)


# Create table of Affect towards Arabs

stargazer(sd_arabs, sd_arabs_right, sd_arabs_left, sd_arabs_center,
          out= "Tables/arab_affect_sd.tex",  style="qje",
          title="Effects of Unity Government on Social Distance towards Arab Citizens of Israel",
          dep.var.labels = c("Social Distance towards Arabs"),
          keep = c("unity", "right_ideology", "left_ideology", 
                   "center_ideology"),
          #omit = c(3:31),
          covariate.labels = c("Unity", "Ideological Right","Ideological Left","Ideological Center",
                               "Unity*Right", "Unity*Left",  "Unity*Center"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
          no.space = T,
          omit.table.layout = "n",
          star.cutoffs = NULL,
          font.size = "small",
          label = "tab:arab_exp_sd",
          star.char = c("", ""),
          add.lines = list(
            c("", "", "", "", ""),
            c("Demographic Controls",  "No", "Yes", "Yes", "Yes")))

# ----------------------------------------------------------------------
# Create Table A26
# ----------------------------------------------------------------------

#Affect towards Haredim (ultra orthodox)

# unity treatment on affect towards haredim
sd_uo <- lm(sdis_uo ~ unity,
            data = main_data)
summary(sd_uo)

sd_uo_control<-lm(sdis_uo ~ unity+ as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                    as.factor(relig) + as.factor(locality),
                  data = main_data)
summary(sd_uo_control)

# unity treatment on affect towards haredim conditional on right
sd_uo_right <- lm(sdis_uo ~ unity*right_ideology+ as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                    as.factor(relig) + as.factor(locality),
                  data = main_data)
summary(sd_uo_right)

# unity treatment on affect towards haredim conditional on left
sd_uo_left <- lm(sdis_uo ~ unity*left_ideology+ as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                   as.factor(relig) + as.factor(locality),
                 data = main_data)
summary(sd_uo_left)

# unity treatment on affect towards haredim conditional on center
sd_uo_center <- lm(sdis_uo ~ unity*center_ideology+ as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                     as.factor(relig) + as.factor(locality),
                   data = main_data)
summary(sd_uo_center)


# Create table of Affect towards Ultra-Orthodox Jews (Haredim)

stargazer(sd_uo, sd_uo_right, sd_uo_left, sd_uo_center,
          out= "Tables/uo_affect_sd.tex",  style="qje",
          title="Effects of Unity Government on Social Distance towards Ultra-Orthodox Jews",
          dep.var.labels = c("Social Distance towards Ultra Orthodox"),
          keep = c("unity", "right_ideology", "left_ideology", 
                   "center_ideology"),
          #omit = c(3:31),
          covariate.labels = c("Unity", "Ideological Right", "Ideological Left",
                               "Ideological Center", "Unity*Right", "Unity*Left",
                               "Unity*Center"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
          no.space = T,
          omit.table.layout = "n",
          star.cutoffs = NULL,
          font.size = "small",
          label = "tab:uo_exp_sd",
          star.char = c("", ""),
          add.lines = list(
            c("", "", "", "", ""),
            c("Demographic Controls",  "No", "Yes", "Yes", "Yes")))


# ----------------------------------------------------------------------
# Create Figure A6
# ----------------------------------------------------------------------
# Sdis with center as left
sd_center <- lm(general_socd_polar ~ unity,
                data = main_data) %>% 
  tidy(., conf.int = T) %>% 
  mutate(.,
         Model = "Social Distance \nCenter as Left")


# Therm with center as left
therm_center <- lm(general_therm_polar ~ unity,
                   data = main_data)%>% 
  tidy(., conf.int = T) %>% 
  mutate(.,
         Model = "Thermometer \nCenter as Left")


# Sdis by vote
sd_by_vote <- lm(prt_sd_polar ~ unity,
                 data = main_data)%>% 
  tidy(., conf.int = T) %>% 
  mutate(.,
         Model = "Social Distance \nBy Vote Choice")

# Therm by vote
therm_by_vote <- lm(prt_therm_polar ~ unity,
                    data = main_data)%>% 
  tidy(., conf.int = T) %>% 
  mutate(.,
         Model = "Thermometer \nBy Vote Choice")


# Sdis 2 party
two_sd <- lm(two_prty_socd_polar ~ unity,
             data = main_data)%>% 
  tidy(., conf.int = T) %>% 
  mutate(.,
         Model = "Social Distance \nTwo Party")

# Therm 2 party
two_therm <- lm(two_prty_therm_polar ~ unity,
                data = main_data)%>% 
  tidy(., conf.int = T) %>% 
  mutate(.,
         Model = "Thermometer \nTwo Party")


# Sdis like in INES
ines_sd <- lm(sd_ines_polar ~ unity,
              data = main_data)%>% 
  tidy(., conf.int = T) %>% 
  mutate(.,
         Model = "Social Distance \nINES Measure")


# Therm like in INES
ines_therm <- lm(therm_ines_polar ~ unity,
                 data = main_data)%>% 
  tidy(., conf.int = T) %>% 
  mutate(.,
         Model = "Thermometer \nINES Measure")

#Socd main parties

main_socd <- lm(main_prty_socd_polar ~ unity,
                data = main_data)%>% 
  tidy(., conf.int = T) %>% 
  mutate(.,
         Model = "Social Distance \nCentral parties")

#Therm main parties

main_therm <- lm(main_prty_therm_polar ~ unity,
                 data = main_data)%>% 
  tidy(., conf.int = T) %>% 
  mutate(.,
         Model = "Thermometer \nCentral parties")

#Socd extreme parties
ext_socd <- lm(ext_prty_socd_polar ~ unity,
               data = main_data)%>% 
  tidy(., conf.int = T) %>% 
  mutate(.,
         Model = "Social Distance \nExtreme parties")

#Therm extreme parties

ext_therm <- lm(ext_prty_therm_polar ~ unity,
                data = main_data)%>% 
  tidy(., conf.int = T) %>% 
  mutate(.,
         Model = "Thermometer \nExtreme parties")

## Generate coefficient plots

sd_rob <- rbind(sd_center, sd_by_vote, main_socd, 
                ext_socd,two_sd,  
                ines_sd) %>% 
  filter(.,
         term == "unity") 

therm_rob <- rbind(therm_center, therm_by_vote, main_therm, 
                   ext_therm, two_therm, 
                   ines_therm) %>% 
  filter(.,
         term == "unity") 

order_sd <- c("Social Distance \nCenter as Left", 
              "Social Distance \nBy Vote Choice",
              "Social Distance \nCentral parties",
              "Social Distance \nExtreme parties",
              "Social Distance \nINES Measure",
              "Social Distance \nTwo Party"
)
ggplot(sd_rob, aes(x = Model, y = estimate)) +
  geom_hline(yintercept = 0, color = "gray50", linetype = 2, size = 0.2) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.4), color = 'darkblue',
                  fill = "white") +
  #coord_flip()+
  labs(x = "",
       y = "Effect Size") +
  scale_x_discrete(limits= order_sd) +
  # caption = "*Coefficients from OLS models, in which outcomes are standardized."
  theme(text = element_text(size = 10, family = "Times"),
        legend.key=element_blank(),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(size = 10),
        plot.caption = element_text(size = 10, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))
ggsave("Figures/sd_exp_rob.pdf", width = 6.5, height = 4)


# ----------------------------------------------------------------------
# Create Figure A7
# ----------------------------------------------------------------------

order_therm <- c("Thermometer \nCenter as Left", 
                 "Thermometer \nBy Vote Choice",
                 "Thermometer \nCentral parties",
                 "Thermometer \nExtreme parties",
                 "Thermometer \nINES Measure",
                 "Thermometer \nTwo Party"
)
ggplot(therm_rob, aes(x = Model, y = estimate)) +
  geom_hline(yintercept = 0, color = "gray50", linetype = 2, size = 0.2) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.4), color = 'darkblue',
                  fill = "white") +
  #coord_flip()+
  labs(x = "",
       y = "Effect Size") +
  scale_x_discrete(limits= order_therm) +
  # caption = "*Coefficients from OLS models, in which outcomes are standardized."
  theme(text = element_text(size = 10, family = "Times"),
        legend.key=element_blank(),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(size = 10),
        plot.caption = element_text(size = 10, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))
ggsave("Figures/therm_exp_rob.pdf", width = 6.5, height = 4)


# ----------------------------------------------------------------------
# Create Table A27
# ----------------------------------------------------------------------

at_soc <- lm(attrition_socd ~ unity,
             data = main_data)
summary(at_soc)

at_therm <- lm(attrition_therm ~ unity,
               data = main_data)
summary(at_therm)

at_soc_cont <- lm(attrition_socd ~ unity + pm_eng +
                    as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                    as.factor(relig) + as.factor(locality) + as.factor(vote),
                  data = main_data)
summary(at_soc_cont)

at_therm_cont <- lm(attrition_therm ~ unity+ pm_eng +
                      as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                      as.factor(relig) + as.factor(locality) + as.factor(vote),
                    data = main_data)
summary(at_therm_cont)

stargazer(at_soc,at_soc_cont, at_therm, at_therm_cont,
          out= "Tables/attrition.tex",  style="qje",
          title="Effect of Treatments on Attrition",
          dep.var.labels = c("Attrition (Social Distance Question)", "Attrition (Thermometer Question)"),
          keep = c("unity", "pm_eng"),
          #omit = c(3:31),
          covariate.labels = c("Unity", "Prime Minister"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
          no.space = T,
          omit.table.layout = "n",
          star.cutoffs = NULL,
          font.size = "small",
          label = "tab:attrition",
          star.char = c("", ""),
          add.lines = list(
            c("", "", "", "", ""),
            c("Demographic Controls",  "No", "Yes", "No", "Yes")))

# ----------------------------------------------------------------------
# Create Table A28
# ----------------------------------------------------------------------
#Plot HTE with SD polarization as outcome

# unity treatment conditional on likud voters
sd_likud <- lm(general_socd_polar_nc ~ unity*likud + as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                 as.factor(relig) + as.factor(locality),
               data = main_data)
summary(sd_likud)

# unity treatment conditional on blue-white voters
sd_bw <- lm(general_socd_polar_nc ~ unity*bw+ as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
              as.factor(relig) + as.factor(locality),
            data = main_data)
summary(sd_bw)

# unity treatment conditional on left
sd_left <- lm(general_socd_polar_nc ~ unity*left_ideology+ as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                as.factor(relig) + as.factor(locality),
              data = main_data)
summary(sd_left)

# unity treatment conditional on right
sd_right <- lm(general_socd_polar_nc ~ unity*right_ideology+ as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                 as.factor(relig) + as.factor(locality),
               data = main_data)
summary(sd_right)

# unity treatment on affect towards right conditional on center 
sd_cen_right <- lm(sdis_right ~ unity*center_ideology+ as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                     as.factor(relig) + as.factor(locality),
                   data = main_data)
summary(sd_cen_right)

# unity treatment on affect towards left conditional on center 
sd_cen_left <- lm(sdis_left ~ unity*center_ideology+ as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                    as.factor(relig) + as.factor(locality),
                  data = main_data)
summary(sd_cen_left)


# Create table of HTE

stargazer(sd_likud, sd_bw, sd_right, sd_left,
          out= "Tables/HTE_experiment_sd.tex",  style="qje",
          title="Effects of Unity Government on Polarization (Social Distance) 
          Conditional on Self-reported vote or Ideology",
          dep.var.labels = c("Polarization (no center)"),
          keep = c("unity", "right_ideology", "left_ideology", 
                   "center_ideology", "bw", "likud"),
          #omit = c(3:31),
          covariate.labels = c("Unity", "Likud", "Blue-White", "Ideological Right", "Ideological Left",
                               "Unity*Likud", "Unity*Blue-White","Unity*Right", "Unity*Left"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
          no.space = T,
          omit.table.layout = "n",
          star.cutoffs = NULL,
          font.size = "small",
          label = "tab:hte_exp_sd",
          star.char = c("", ""),
          add.lines = list(
            c("", "", "", "", ""),
            c("Demographic Controls",  "Yes", "Yes", "Yes", "Yes")))


# ----------------------------------------------------------------------
# Create Table A29
# ----------------------------------------------------------------------

# unity treatment conditional on likud voters
sd_netan <- lm(general_socd_polar_nc ~ unity*pm_eng,
               data = main_data)
summary(sd_netan)

sd_netan_cov <- lm(general_socd_polar_nc ~ unity*pm_eng + as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                     as.factor(relig) + as.factor(locality),
                   data = main_data)
summary(sd_netan_cov)


# Create table of HTE

stargazer(sd_netan, sd_netan_cov, 
          out= "Tables/HTE_netan_experiment_sd.tex",  style="qje",
          title="Effects of Unity Government on Polarization (Social Distance) 
          Conditional on PM Treatment",
          dep.var.labels = c("Polarization (no center)"),
          keep = c("unity", "pm_engNetanyahu", "unity:pm_engNetanyahu"),
          #omit = c(3:31),
          covariate.labels = c("Unity", "Netanyahu", "Unity*Netanyahu"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
          no.space = T,
          omit.table.layout = "n",
          star.cutoffs = NULL,
          font.size = "small",
          label = "tab:hte_netan_exp_sd",
          star.char = c("", ""),
          add.lines = list(
            c("", "", "", "", ""),
            c("Demographic Controls",  "No", "Yes")))


# ----------------------------------------------------------------------
# Create Table A30
# ----------------------------------------------------------------------

# Main model using SD polarization without center

# Polar base
netanyahu_1 <- lm(general_socd_polar_nc ~ pm_eng,
                  data = main_data)
summary(netanyahu_1)

# Polar base + demographics
netanyahu_2 <- lm(general_socd_polar_nc ~ pm_eng + left_right_1 +
                    as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                    as.factor(relig) + as.factor(locality) + as.factor(vote),
                  data = main_data)
summary(netanyahu_2)

# Polar base interaction of netanyahu*left_wing + demographics
netanyahu_3 <- lm(general_socd_polar_nc ~ pm_eng*left_ideology +
                    as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                    as.factor(relig) + as.factor(locality) + as.factor(vote),
                  data = main_data)
summary(netanyahu_3)

# Polar base interaction of netanyahu*right_wing + demographics
netanyahu_4 <- lm(general_socd_polar_nc ~ pm_eng*right_ideology +
                    as.factor(sex) + as.factor(age) + as.factor(ethnicity) +
                    as.factor(relig) + as.factor(locality) + as.factor(vote),
                  data = main_data)
summary(netanyahu_4)


# Create table of main results
stargazer(netanyahu_1, netanyahu_2, netanyahu_3, netanyahu_4,
          out= "Tables/netanyahu_sd.tex",  style="qje",
          title="Effects of Netanyahu Treatment on Polarization (Social Distance)",
          dep.var.labels = "Polarization",
          
          
          keep = c("pm_eng", "left_ideology","right_ideology", 
                   "pm_eng*left_ideology", 
                   "pm_eng*right_ideology"),
          #omit = c(3:31),
          covariate.labels = c("PM", "Left","Right",
                               "PM*Left",  "PM*Right"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
          no.space = T,
          omit.table.layout = "n",
          star.cutoffs = NULL,
          font.size = "small",
          label = "tab:netanyahu_sd",
          star.char = c("", "", ""),
          add.lines = list(
            c("", "", "", "", "", "" ),
            c("Demographic Controls",  "No", "Yes", "Yes", "Yes"),
            c("Center Voters", "No", "No", "No", "No")))

